%% Figure C1: Labor Market Earnings: Data and Estimated Lognormal Distributions
% Code for generating histogram and density of Earnings by education and 
% occupation 

clear all

filename = 'Inc_educ_routine_cognitive_5';
sheet = 'Sheet1';

mu=6:0.001:8.5;
sigma=0.6:0.001:1.1;
l_mu=length(mu);
l_sigma=length(sigma);

x = xlsread(filename,sheet,'A2:A51'); %ln income
y1 = xlsread(filename,sheet,'B2:B51'); %education 0 routine
y2 = xlsread(filename,sheet,'C2:C51'); %education 1-4 routine
y3 = xlsread(filename,sheet,'D2:D51'); %education 5-8 routine
y4 = xlsread(filename,sheet,'E2:E51'); %education 9-11 routine
y5 = xlsread(filename,sheet,'F2:F51'); %education 12+ cognitive

y6 = xlsread(filename,sheet,'G2:G51'); %education 0 cognitive
y7 = xlsread(filename,sheet,'H2:H51'); %education 1-4 cognitive
y8 = xlsread(filename,sheet,'I2:I51'); %education 5-8 cognitive
y9 = xlsread(filename,sheet,'J2:J51'); %education 9-11 cognitive
y10 = xlsread(filename,sheet,'K2:K51'); %education 12+ cognitive

for i=1:l_mu
    for j=1:l_sigma
        z(i,j,:)=normpdf(x,mu(i),sigma(j));
        y1_matrix(i,j,:)=transpose(y1);
        dif_sum1(i,j)=sum(abs(z(i,j,:)-y1_matrix(i,j,:)));
        y2_matrix(i,j,:)=transpose(y2);
        dif_sum2(i,j)=sum(abs(z(i,j,:)-y2_matrix(i,j,:)));
        y3_matrix(i,j,:)=transpose(y3);
        dif_sum3(i,j)=sum(abs(z(i,j,:)-y3_matrix(i,j,:)));
        y4_matrix(i,j,:)=transpose(y4);
        dif_sum4(i,j)=sum(abs(z(i,j,:)-y4_matrix(i,j,:)));
        y5_matrix(i,j,:)=transpose(y5);
        dif_sum5(i,j)=sum(abs(z(i,j,:)-y5_matrix(i,j,:)));
        y6_matrix(i,j,:)=transpose(y6);
        dif_sum6(i,j)=sum(abs(z(i,j,:)-y6_matrix(i,j,:)));
        y7_matrix(i,j,:)=transpose(y7);
        dif_sum7(i,j)=sum(abs(z(i,j,:)-y7_matrix(i,j,:)));
        y8_matrix(i,j,:)=transpose(y8);
        dif_sum8(i,j)=sum(abs(z(i,j,:)-y8_matrix(i,j,:)));
        y9_matrix(i,j,:)=transpose(y9);
        dif_sum9(i,j)=sum(abs(z(i,j,:)-y9_matrix(i,j,:)));
        y10_matrix(i,j,:)=transpose(y10);
        dif_sum10(i,j)=sum(abs(z(i,j,:)-y10_matrix(i,j,:)));
    end
end


[min_val1,idx1]=min(dif_sum1(:));
[row1,col1]=ind2sub(size(dif_sum1),idx1);
min_mu_1 = mu(row1);
min_sigma_1 = sigma(col1);

[min_val2,idx2]=min(dif_sum2(:));
[row2,col2]=ind2sub(size(dif_sum2),idx2);
min_mu_2 = mu(row2);
min_sigma_2 = sigma(col2);

[min_val3,idx3]=min(dif_sum3(:));
[row3,col3]=ind2sub(size(dif_sum3),idx3);
min_mu_3 = mu(row3);
min_sigma_3 = sigma(col3);

[min_val4,idx4]=min(dif_sum4(:));
[row4,col4]=ind2sub(size(dif_sum4),idx4);
min_mu_4 = mu(row4);
min_sigma_4 = sigma(col4);

[min_val5,idx5]=min(dif_sum5(:));
[row5,col5]=ind2sub(size(dif_sum5),idx5);
min_mu_5 = mu(row5);
min_sigma_5 = sigma(col5);

[min_val6,idx6]=min(dif_sum6(:));
[row6,col6]=ind2sub(size(dif_sum6),idx6);
min_mu_6 = mu(row6);
min_sigma_6 = sigma(col6);

[min_val7,idx7]=min(dif_sum7(:));
[row7,col7]=ind2sub(size(dif_sum7),idx7);
min_mu_7 = mu(row7);
min_sigma_7 = sigma(col7);

[min_val8,idx8]=min(dif_sum8(:));
[row8,col8]=ind2sub(size(dif_sum8),idx8);
min_mu_8 = mu(row8);
min_sigma_8 = sigma(col8);

[min_val9,idx9]=min(dif_sum9(:));
[row9,col9]=ind2sub(size(dif_sum9),idx9);
min_mu_9 = mu(row9);
min_sigma_9 = sigma(col9);

[min_val10,idx10]=min(dif_sum10(:));
[row10,col10]=ind2sub(size(dif_sum10),idx10);
min_mu_10 = mu(row10);
min_sigma_10 = sigma(col10);

z1 = normpdf(x,min_mu_1,min_sigma_1);
z2 = normpdf(x,min_mu_2,min_sigma_2);
z3 = normpdf(x,min_mu_3,min_sigma_3);
z4 = normpdf(x,min_mu_4,min_sigma_4); 
z5 = normpdf(x,min_mu_5,min_sigma_5);
z6 = normpdf(x,min_mu_6,min_sigma_6);
z7 = normpdf(x,min_mu_7,min_sigma_7);
z8 = normpdf(x,min_mu_8,min_sigma_8);
z9 = normpdf(x,min_mu_9,min_sigma_9);
z10 = normpdf(x,min_mu_10,min_sigma_10);

figure (1)
plot(x,y1,'Color',[0 0 0],'LineWidth', 2.5)
hold on
plot(x,z1,'Color',[0 0 0],'LineStyle','--','LineWidth', 2.5)
hold on 
plot(x,y2,'Color',[0.3 0.3 0.3],'LineWidth', 2.5)
hold on
plot(x,z2,'Color',[0.3 0.3 0.3],'LineStyle','--','LineWidth', 2.5)
hold on 
plot(x,y3,'Color',[0.5 0.5 0.5],'LineWidth', 2.5)
hold on
plot(x,z3,'Color',[0.5 0.5 0.5],'LineStyle','--','LineWidth', 2.5)
hold on 
plot(x,y4,'Color',[0.7 0.7 0.7],'LineWidth', 2.5)
hold on
plot(x,z4,'Color',[0.7 0.7 0.7],'LineStyle','--','LineWidth', 2.5)
hold on 
plot(x,y5,'Color',[0.85 0.85 0.85],'LineWidth', 2.5)
hold on
plot(x,z5,'Color',[0.85 0.85 0.85],'LineStyle','--','LineWidth', 2.5)
xlabel('Ln Income') 
ylabel('Frequency')
legend({['School=0'],['Normal' 10 '(mu=6.22 & sigma=0.67)'],['School=1-4'],['Normal' 10 '(mu=6.67 & sigma=0.72)'],['School=5-8'],['Normal' 10 '(mu=7.01 & sigma=0.78)'],['School=9-11'],['Normal' 10 '(mu=7.34 & sigma=0.82)'],['School=12+'],['Normal' 10 '(mu=7.83 & sigma=0.83)']},'FontWeight','bold', 'Location','best','Orientation','vertical')

figure (2)
plot(x,y6,'Color',[0 0 0],'LineWidth', 2.5)
hold on
plot(x,z6,'Color',[0 0 0],'LineStyle','--','LineWidth', 2.5)
hold on 
plot(x,y7,'Color',[0.3 0.3 0.3],'LineWidth', 2.5)
hold on
plot(x,z7,'Color',[0.3 0.3 0.3],'LineStyle','--','LineWidth', 2.5)
hold on 
plot(x,y8,'Color',[0.5 0.5 0.5],'LineWidth', 2.5)
hold on
plot(x,z8,'Color',[0.5 0.5 0.5],'LineStyle','--','LineWidth', 2.5)
hold on 
plot(x,y9,'Color',[0.7 0.7 0.7],'LineWidth', 2.5)
hold on
plot(x,z9,'Color',[0.7 0.7 0.7],'LineStyle','--','LineWidth', 2.5)
hold on 
plot(x,y10,'Color',[0.85 0.85 0.85],'LineWidth', 2.5)
hold on
plot(x,z10,'Color',[0.85 0.85 0.85],'LineStyle','--','LineWidth', 2.5)
ylim([0 0.5])
xlabel('Ln Income') 
ylabel('Frequency')
legend({['School=0'],['Normal' 10 '(mu=6.63 & sigma=1.00)'],['School=1-4'],['Normal' 10 '(mu=7.39 & sigma=1.03)'],['School=5-8'],['Normal' 10 '(mu=7.59 & sigma=0.93)'],['School=9-11'],['Normal' 10 '(mu=7.46 & sigma=1.04)'],['School=12+'],['Normal' 10 '(mu=8.29 & sigma=0.95)']},'FontWeight','bold', 'Location','best','Orientation','vertical')
